#### Meta Data ####

# article: "The Ambivalent Effect of Autocratization on Domestic Terrorism" 
# journal: "Studies in Conflict & Terrorism"
# authors: "Lott, Lars; Croissant, Aurel; Trinn. Christoph"
# date: 2023-10-02
# written under "R version 4.3.1 (2023-06-16 ucrt)"


#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# set working directory
getwd()

#### Preliminaries ####

library(MASS)
library(multiwayvcov)
library(ggeffects)
library(tidyverse)
library(ggpubr)
library(here)
library(haven)
library(imputeTS)
library(modelsummary)
library(data.table)



#### functions ####

panel_match_path <- "PanelMatch/R/"

source(file.path(panel_match_path, "data.R"))
source(file.path(panel_match_path, "DisplayTreatment.R"))
source(file.path(panel_match_path, "exact_match.R"))
source(file.path(panel_match_path, "listwise_delete.R"))
source(file.path(panel_match_path, "matched_set_obj.R"))
source(file.path(panel_match_path, "matched_set_R.R"))
source(file.path(panel_match_path, "network_and_caliper.R"))
source(file.path(panel_match_path, "PanelEstimate.R"))
source(file.path(panel_match_path, "PanelEstimateObject.R"))
source(file.path(panel_match_path, "PanelMatch.R"))
source(file.path(panel_match_path, "PanelMatch-Package.R"))
source(file.path(panel_match_path, "PE_helpers.R"))
source(file.path(panel_match_path, "PE_lower_level.R"))
source(file.path(panel_match_path, "pm_helpers_r.R"))
source(file.path(panel_match_path, "RcppExports.R"))
source(file.path(panel_match_path, "utilities.R"))



#### Load Dataset ####

vdem <- readRDS("data/vdem12.rds") 

summary(vdem$aut_ep_outcome_agg)

#### Define Treatment ####

vdem <- vdem %>%
  mutate(autocratization_treatment = aut_ep,
         democratic_breakdown = ifelse(aut_ep_outcome_agg == 1, 1, 0), 
         democratic_recession = ifelse(aut_ep_outcome_agg == 2, 1, 0), 
         regressed_autocracy = ifelse(aut_ep_outcome_agg == 3, 1, 0))

table(vdem$autocratization_treatment)
table(vdem$democratic_breakdown)
table(vdem$democratic_recession)
table(vdem$regressed_autocracy)

vdem <- vdem %>%
  mutate(democratic_breakdown_recession = ifelse(democratic_breakdown == 1 |democratic_recession == 1, 1,0 ))

table(vdem$democratic_breakdown_recession)



vdem <- vdem %>%
  mutate(log_count_attacks_dom = log(count_attacks_dom +0.01),
         log_count_attacks_int = log(count_attacks_int +0.01), 
         log_count_attacks = log(count_attacks +0.01),
         log_sum_nkill_dom = log(sum_nkill_dom +0.01),
         log_sum_nkill_int = log(sum_nkill_int +0.01),
         log_sum_nkill = log(sum_nkill +0.01)) %>%
  filter(year >= 1966)


vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(log_count_attacks_dom = na_interpolation(log_count_attacks_dom, option = "spline",  maxgap = 1), 
         log_sum_nkill_dom = na_interpolation(log_sum_nkill_dom, option = "spline",  maxgap = 1), 
         count_attacks_dom = na_interpolation(count_attacks_dom, option = "spline",  maxgap = 1), 
         sum_nkill_dom = na_interpolation(sum_nkill_dom, option = "spline",  maxgap = 1), 
         log_count_attacks_dom = ifelse(log_count_attacks_dom<0, 0, log_count_attacks_dom),
         log_sum_nkill_dom = ifelse(log_sum_nkill_dom<0, 0, log_sum_nkill_dom),
         count_attacks_dom = ifelse(count_attacks_dom<0, 0, count_attacks_dom),
         sum_nkill_dom = ifelse(sum_nkill_dom<0, 0, sum_nkill_dom))

vdem <- vdem %>%
  dplyr::select(country_id, country_name, year, aut_ep, 
                aut_ep_id,  log_count_attacks_dom, log_count_attacks_int, log_count_attacks, 
                log_sum_nkill_dom, log_sum_nkill_int, log_sum_nkill, count_attacks_dom, count_attacks_int, count_attacks, 
                sum_nkill_dom, sum_nkill_int, sum_nkill,e_gdppc, e_pop,
                v2x_polyarchy, milper_extrapol, eprratio_extrapol, democratic_recession, democratic_breakdown, 
                civil_war_ucdp, v2x_regime, aut_ep_outcome, aut_ep_outcome_agg, 
                v2caviol, dem_ep, dem_ep_id, dem_ep_outcome_agg, dem_ep_outcome, e_regionpol, e_regionpol_6C, 
                autocratization_treatment, regressed_autocracy, democratic_breakdown_recession, ma_count_attacks_dom, ma_sum_nkill_dom, democracy_stock, 
                v2x_veracc, v2x_diagacc, v2x_horacc, e_pelifeex, e_peaveduc, e_area) %>%
  filter(year <=2020) %>%
  group_by(country_id) %>%
  fill(e_area)

vdem <- vdem %>%
  drop_na(autocratization_treatment) %>%
  mutate(e_gdppc_sqaured = e_gdppc^2, 
         e_area_log = log10(e_area))

summary(vdem$e_area_log)
summary(vdem$e_gdppc_sqaured)

#------------------------------------------------------------------------------#
#----------------------------- Figures G3 and G6 ------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)
vdem <- as.data.frame(vdem)
class(vdem)


# Define Lags and leads for all

lag_2 <- c(1:2)
lead_5 <- c(0:5)

vdem <- vdem %>%
  dplyr::select(-c(country_name, aut_ep_id))

vdem$country_id <- as.integer(vdem$country_id)

#### No refinement #### 

PM_M1 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_treatment", refinement.method = "none", 
                        data = vdem, match.missing = FALSE, 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE, placebo.test = T)

M1 <- PanelEstimate(sets = PM_M1, data = vdem)

M1_ci09 <- PanelEstimate(sets = PM_M1, data = vdem, confidence.level = .9)

M1_placebo <- placebo_test(PM_M1,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M1_placebo_t_2_est <- M1_placebo$estimates[1]
M1_placebo_t_2_high <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.975)
M1_placebo_t_2_low <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.025)

M1_placebo_t_2_high_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.95)
M1_placebo_t_2_low_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.05)


M1_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M1_placebo_t_2_est,
                              conf_low_95 = M1_placebo_t_2_low,
                              conf_low_90 = M1_placebo_t_2_low_90,
                              conf_high_95 = M1_placebo_t_2_high, 
                              conf_high_90 = M1_placebo_t_2_high_90)

#95% CIs ##

M1_pred <- summary(M1)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M1_pred_ci09 <- summary(M1_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_1 <- as.data.frame(M1_pred)
data_figure_1_09 <- as.data.frame(M1_pred_ci09)

data_figure_1 <- cbind(data_figure_1, data_figure_1_09[,3:4])

data_figure_1 <- data_figure_1  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_1 <- rbind(data_figure_1, M1_placebo_data)

data_figure_1[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Autocraitzation", 
       y = "Estimated effect of Autocratization", 
       title = "No Covariates: Terrorist Attacks")

#### Mahalanobis #### 

PM_M2 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "autocratization_treatment", refinement.method = "mahalanobis", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M2 <- PanelEstimate(sets = PM_M2, data = vdem)

M2_ci09 <- PanelEstimate(sets = PM_M2, data = vdem, confidence.level = .9)

M2_placebo <- placebo_test(PM_M2,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M2_placebo_t_2_est <- M2_placebo$estimates[1]
M2_placebo_t_2_high <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.975)
M2_placebo_t_2_low <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.025)

M2_placebo_t_2_high_90 <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.95)
M2_placebo_t_2_low_90 <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.05)


M2_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M2_placebo_t_2_est,
                              conf_low_95 = M2_placebo_t_2_low,
                              conf_low_90 = M2_placebo_t_2_low_90,
                              conf_high_95 = M2_placebo_t_2_high, 
                              conf_high_90 = M2_placebo_t_2_high_90)

#95% CIs ##

M2_pred <- summary(M2)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M2_pred_ci09 <- summary(M2_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_2 <- as.data.frame(M2_pred)
data_figure_2_09 <- as.data.frame(M2_pred_ci09)

data_figure_2 <- cbind(data_figure_2, data_figure_2_09[,3:4])

data_figure_2 <- data_figure_2  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_2 <- rbind(data_figure_2, M2_placebo_data)

data_figure_2[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1b <- ggplot(data = data_figure_2, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Autocraitzation", 
       y = "Estimated effect of Autocratization", 
       title = "Mahalanobis: Terrorist Attacks")


#### PS Match #### 

PM_M3 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "autocratization_treatment", refinement.method = "ps.match", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M3 <- PanelEstimate(sets = PM_M3, data = vdem)

M3_ci09 <- PanelEstimate(sets = PM_M3, data = vdem, confidence.level = .9)

M3_placebo <- placebo_test(PM_M3,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M3_placebo_t_2_est <- M3_placebo$estimates[1]
M3_placebo_t_2_high <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.975)
M3_placebo_t_2_low <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.025)

M3_placebo_t_2_high_90 <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.95)
M3_placebo_t_2_low_90 <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.05)


M3_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M3_placebo_t_2_est,
                              conf_low_95 = M3_placebo_t_2_low,
                              conf_low_90 = M3_placebo_t_2_low_90,
                              conf_high_95 = M3_placebo_t_2_high, 
                              conf_high_90 = M3_placebo_t_2_high_90)

#95% CIs ##

M3_pred <- summary(M3)

M3_pred <- M3_pred$summary

M3_pred <- cbind(M3_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M3_pred_ci09 <- summary(M3_ci09)

M3_pred_ci09 <- M3_pred_ci09$summary

M3_pred_ci09 <- cbind(M3_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_3 <- as.data.frame(M3_pred)
data_figure_3_09 <- as.data.frame(M3_pred_ci09)

data_figure_3 <- cbind(data_figure_3, data_figure_3_09[,3:4])

data_figure_3 <- data_figure_3  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_3 <- rbind(data_figure_3, M3_placebo_data)

data_figure_3[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1c <- ggplot(data = data_figure_3, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Autocraitzation", 
       y = "Estimated effect of Autocratization", 
       title = "PS Match: Terrorist Attacks")




#### CBPS Match #### 

PM_M4 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "autocratization_treatment", refinement.method = "CBPS.match", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M4 <- PanelEstimate(sets = PM_M4, data = vdem)

M4_ci09 <- PanelEstimate(sets = PM_M4, data = vdem, confidence.level = .9)

M4_placebo <- placebo_test(PM_M4,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M4_placebo_t_2_est <- M4_placebo$estimates[1]
M4_placebo_t_2_high <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.975)
M4_placebo_t_2_low <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.025)

M4_placebo_t_2_high_90 <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.95)
M4_placebo_t_2_low_90 <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.05)


M4_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M4_placebo_t_2_est,
                              conf_low_95 = M4_placebo_t_2_low,
                              conf_low_90 = M4_placebo_t_2_low_90,
                              conf_high_95 = M4_placebo_t_2_high, 
                              conf_high_90 = M4_placebo_t_2_high_90)

#95% CIs ##

M4_pred <- summary(M4)

M4_pred <- M4_pred$summary

M4_pred <- cbind(M4_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M4_pred_ci09 <- summary(M4_ci09)

M4_pred_ci09 <- M4_pred_ci09$summary

M4_pred_ci09 <- cbind(M4_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_4 <- as.data.frame(M4_pred)
data_figure_4_09 <- as.data.frame(M4_pred_ci09)

data_figure_4 <- cbind(data_figure_4, data_figure_4_09[,3:4])

data_figure_4 <- data_figure_4  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_4 <- rbind(data_figure_4, M4_placebo_data)

data_figure_4[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1d <- ggplot(data = data_figure_4, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Autocraitzation", 
       y = "Estimated effect of Autocratization", 
       title = " CBPS Match: Terrorist Attacks")



#### CBPS Weight #### 

PM_M5 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "autocratization_treatment", refinement.method = "CBPS.weight", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M5<- PanelEstimate(sets = PM_M5, data = vdem)

M5_ci09 <- PanelEstimate(sets = PM_M5, data = vdem, confidence.level = .9)

M5_placebo <- placebo_test(PM_M5,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M5_placebo_t_2_est <- M5_placebo$estimates[1]
M5_placebo_t_2_high <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.975)
M5_placebo_t_2_low <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.025)

M5_placebo_t_2_high_90 <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.95)
M5_placebo_t_2_low_90 <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.05)


M5_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M5_placebo_t_2_est,
                              conf_low_95 = M5_placebo_t_2_low,
                              conf_low_90 = M5_placebo_t_2_low_90,
                              conf_high_95 = M5_placebo_t_2_high, 
                              conf_high_90 = M5_placebo_t_2_high_90)

#95% CIs ##

M5_pred <- summary(M5)

M5_pred <- M5_pred$summary

M5_pred <- cbind(M5_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M5_pred_ci09 <- summary(M5_ci09)

M5_pred_ci09 <- M5_pred_ci09$summary

M5_pred_ci09 <- cbind(M5_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_5 <- as.data.frame(M5_pred)
data_figure_5_09 <- as.data.frame(M5_pred_ci09)

data_figure_5 <- cbind(data_figure_5, data_figure_5_09[,3:4])

data_figure_5 <- data_figure_5  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_5 <- rbind(data_figure_5, M5_placebo_data)

data_figure_5[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1e <- ggplot(data = data_figure_5, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Autocraitzation", 
       y = "Estimated effect of Autocratization", 
       title = "CBPS Weight: Terrorist Attacks")

##### GG Arrange Plots #####

ggarrange(f1a, f1b, f1c, f1d, f1e, nrow = 5)

ggsave("outputs/Figure_G6.png", units= c("cm"), width = 20, height = 35, dpi = 900)
ggsave("outputs/Figure_G6.pdf", units= c("cm"), width = 20, height = 35, dpi = 900)


#------------------------------------------------------------------------------#
#------------------------------ Figure G4  -------------------------------------#
#------------------------------------------------------------------------------#

summary(vdem$democratic_breakdown_recession)
table(vdem$democratic_breakdown_recession)


#### No refinement #### 

PM_M1 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "democratic_breakdown_recession", refinement.method = "none", 
                    data = vdem, match.missing = FALSE, 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M1 <- PanelEstimate(sets = PM_M1, data = vdem)

M1_ci09 <- PanelEstimate(sets = PM_M1, data = vdem, confidence.level = .9)

M1_placebo <- placebo_test(PM_M1,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M1_placebo_t_2_est <- M1_placebo$estimates[1]
M1_placebo_t_2_high <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.975)
M1_placebo_t_2_low <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.025)

M1_placebo_t_2_high_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.95)
M1_placebo_t_2_low_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.05)


M1_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M1_placebo_t_2_est,
                              conf_low_95 = M1_placebo_t_2_low,
                              conf_low_90 = M1_placebo_t_2_low_90,
                              conf_high_95 = M1_placebo_t_2_high, 
                              conf_high_90 = M1_placebo_t_2_high_90)

#95% CIs ##

M1_pred <- summary(M1)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M1_pred_ci09 <- summary(M1_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_1 <- as.data.frame(M1_pred)
data_figure_1_09 <- as.data.frame(M1_pred_ci09)

data_figure_1 <- cbind(data_figure_1, data_figure_1_09[,3:4])

data_figure_1 <- data_figure_1  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_1 <- rbind(data_figure_1, M1_placebo_data)

data_figure_1[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Democratic Breakdown/Recession", 
       y = "Estimated effect of Democratic Breakdown/Recession", 
       title = "No Covariates: Terrorist Attacks")

#### Mahalanobis #### 

PM_M2 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "democratic_breakdown_recession", refinement.method = "mahalanobis", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M2 <- PanelEstimate(sets = PM_M2, data = vdem)

M2_ci09 <- PanelEstimate(sets = PM_M2, data = vdem, confidence.level = .9)

M2_placebo <- placebo_test(PM_M2,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M2_placebo_t_2_est <- M2_placebo$estimates[1]
M2_placebo_t_2_high <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.975)
M2_placebo_t_2_low <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.025)

M2_placebo_t_2_high_90 <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.95)
M2_placebo_t_2_low_90 <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.05)


M2_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M2_placebo_t_2_est,
                              conf_low_95 = M2_placebo_t_2_low,
                              conf_low_90 = M2_placebo_t_2_low_90,
                              conf_high_95 = M2_placebo_t_2_high, 
                              conf_high_90 = M2_placebo_t_2_high_90)

#95% CIs ##

M2_pred <- summary(M2)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M2_pred_ci09 <- summary(M2_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_2 <- as.data.frame(M2_pred)
data_figure_2_09 <- as.data.frame(M2_pred_ci09)

data_figure_2 <- cbind(data_figure_2, data_figure_2_09[,3:4])

data_figure_2 <- data_figure_2  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_2 <- rbind(data_figure_2, M2_placebo_data)

data_figure_2[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1b <- ggplot(data = data_figure_2, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Democratic Breakdown/Recession", 
       y = "Estimated effect of Democratic Breakdown/Recession", 
       title = "Mahalanobis: Terrorist Attacks")


#### PS Match #### 

PM_M3 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "democratic_breakdown_recession", refinement.method = "ps.match", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M3 <- PanelEstimate(sets = PM_M3, data = vdem)

M3_ci09 <- PanelEstimate(sets = PM_M3, data = vdem, confidence.level = .9)

M3_placebo <- placebo_test(PM_M3,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M3_placebo_t_2_est <- M3_placebo$estimates[1]
M3_placebo_t_2_high <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.975)
M3_placebo_t_2_low <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.025)

M3_placebo_t_2_high_90 <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.95)
M3_placebo_t_2_low_90 <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.05)


M3_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M3_placebo_t_2_est,
                              conf_low_95 = M3_placebo_t_2_low,
                              conf_low_90 = M3_placebo_t_2_low_90,
                              conf_high_95 = M3_placebo_t_2_high, 
                              conf_high_90 = M3_placebo_t_2_high_90)

#95% CIs ##

M3_pred <- summary(M3)

M3_pred <- M3_pred$summary

M3_pred <- cbind(M3_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M3_pred_ci09 <- summary(M3_ci09)

M3_pred_ci09 <- M3_pred_ci09$summary

M3_pred_ci09 <- cbind(M3_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_3 <- as.data.frame(M3_pred)
data_figure_3_09 <- as.data.frame(M3_pred_ci09)

data_figure_3 <- cbind(data_figure_3, data_figure_3_09[,3:4])

data_figure_3 <- data_figure_3  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_3 <- rbind(data_figure_3, M3_placebo_data)

data_figure_3[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1c <- ggplot(data = data_figure_3, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Democratic Breakdown/Recession", 
       y = "Estimated effect of Democratic Breakdown/Recession", 
       title = "PS Match: Terrorist Attacks")




#### CBPS Match #### 

PM_M4 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "democratic_breakdown_recession", refinement.method = "CBPS.match", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M4 <- PanelEstimate(sets = PM_M4, data = vdem)

M4_ci09 <- PanelEstimate(sets = PM_M4, data = vdem, confidence.level = .9)

M4_placebo <- placebo_test(PM_M4,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M4_placebo_t_2_est <- M4_placebo$estimates[1]
M4_placebo_t_2_high <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.975)
M4_placebo_t_2_low <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.025)

M4_placebo_t_2_high_90 <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.95)
M4_placebo_t_2_low_90 <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.05)


M4_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M4_placebo_t_2_est,
                              conf_low_95 = M4_placebo_t_2_low,
                              conf_low_90 = M4_placebo_t_2_low_90,
                              conf_high_95 = M4_placebo_t_2_high, 
                              conf_high_90 = M4_placebo_t_2_high_90)

#95% CIs ##

M4_pred <- summary(M4)

M4_pred <- M4_pred$summary

M4_pred <- cbind(M4_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M4_pred_ci09 <- summary(M4_ci09)

M4_pred_ci09 <- M4_pred_ci09$summary

M4_pred_ci09 <- cbind(M4_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_4 <- as.data.frame(M4_pred)
data_figure_4_09 <- as.data.frame(M4_pred_ci09)

data_figure_4 <- cbind(data_figure_4, data_figure_4_09[,3:4])

data_figure_4 <- data_figure_4  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_4 <- rbind(data_figure_4, M4_placebo_data)

data_figure_4[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1d <- ggplot(data = data_figure_4, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Democratic Breakdown/Recession", 
       y = "Estimated effect of Democratic Breakdown/Recession", 
       title = " CBPS Match: Terrorist Attacks")



#### CBPS Weight #### 

PM_M5 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "democratic_breakdown_recession", refinement.method = "CBPS.weight", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M5<- PanelEstimate(sets = PM_M5, data = vdem)

M5_ci09 <- PanelEstimate(sets = PM_M5, data = vdem, confidence.level = .9)

M5_placebo <- placebo_test(PM_M5,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M5_placebo_t_2_est <- M5_placebo$estimates[1]
M5_placebo_t_2_high <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.975)
M5_placebo_t_2_low <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.025)

M5_placebo_t_2_high_90 <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.95)
M5_placebo_t_2_low_90 <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.05)


M5_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M5_placebo_t_2_est,
                              conf_low_95 = M5_placebo_t_2_low,
                              conf_low_90 = M5_placebo_t_2_low_90,
                              conf_high_95 = M5_placebo_t_2_high, 
                              conf_high_90 = M5_placebo_t_2_high_90)

#95% CIs ##

M5_pred <- summary(M5)

M5_pred <- M5_pred$summary

M5_pred <- cbind(M5_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M5_pred_ci09 <- summary(M5_ci09)

M5_pred_ci09 <- M5_pred_ci09$summary

M5_pred_ci09 <- cbind(M5_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_5 <- as.data.frame(M5_pred)
data_figure_5_09 <- as.data.frame(M5_pred_ci09)

data_figure_5 <- cbind(data_figure_5, data_figure_5_09[,3:4])

data_figure_5 <- data_figure_5  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_5 <- rbind(data_figure_5, M5_placebo_data)

data_figure_5[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1e <- ggplot(data = data_figure_5, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Democratic Breakdown/Recession", 
       y = "Estimated effect of Democratic Breakdown/Recession", 
       title = "CBPS Weight: Terrorist Attacks")

##### GG Arrange Plots #####

ggarrange(f1a, f1b, f1c, f1d, f1e, nrow = 5)

ggsave("outputs/Figure_G4.png", units= c("cm"), width = 20, height = 35, dpi = 900)
ggsave("outputs/Figure_G4.pdf", units= c("cm"), width = 20, height = 35, dpi = 900)


#------------------------------------------------------------------------------#
#------------------------------ Figure_G5    ----------------------------------#
#------------------------------------------------------------------------------#

summary(vdem$regressed_autocracy)
table(vdem$regressed_autocracy)


#### No refinement #### 

PM_M1 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "regressed_autocracy", refinement.method = "none", 
                    data = vdem, match.missing = FALSE, 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M1 <- PanelEstimate(sets = PM_M1, data = vdem)

M1_ci09 <- PanelEstimate(sets = PM_M1, data = vdem, confidence.level = .9)

M1_placebo <- placebo_test(PM_M1,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M1_placebo_t_2_est <- M1_placebo$estimates[1]
M1_placebo_t_2_high <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.975)
M1_placebo_t_2_low <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.025)

M1_placebo_t_2_high_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.95)
M1_placebo_t_2_low_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.05)


M1_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M1_placebo_t_2_est,
                              conf_low_95 = M1_placebo_t_2_low,
                              conf_low_90 = M1_placebo_t_2_low_90,
                              conf_high_95 = M1_placebo_t_2_high, 
                              conf_high_90 = M1_placebo_t_2_high_90)

#95% CIs ##

M1_pred <- summary(M1)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M1_pred_ci09 <- summary(M1_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_1 <- as.data.frame(M1_pred)
data_figure_1_09 <- as.data.frame(M1_pred_ci09)

data_figure_1 <- cbind(data_figure_1, data_figure_1_09[,3:4])

data_figure_1 <- data_figure_1  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_1 <- rbind(data_figure_1, M1_placebo_data)

data_figure_1[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Autocratic Consolidation", 
       y = "Estimated effect of Autocratic Consolidation", 
       title = "No Covariates: Terrorist Attacks")

#### Mahalanobis #### 

PM_M2 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "regressed_autocracy", refinement.method = "mahalanobis", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M2 <- PanelEstimate(sets = PM_M2, data = vdem)

M2_ci09 <- PanelEstimate(sets = PM_M2, data = vdem, confidence.level = .9)

M2_placebo <- placebo_test(PM_M2,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M2_placebo_t_2_est <- M2_placebo$estimates[1]
M2_placebo_t_2_high <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.975)
M2_placebo_t_2_low <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.025)

M2_placebo_t_2_high_90 <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.95)
M2_placebo_t_2_low_90 <- quantile(M2_placebo$bootstrapped.estimates[ ,1],0.05)


M2_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M2_placebo_t_2_est,
                              conf_low_95 = M2_placebo_t_2_low,
                              conf_low_90 = M2_placebo_t_2_low_90,
                              conf_high_95 = M2_placebo_t_2_high, 
                              conf_high_90 = M2_placebo_t_2_high_90)

#95% CIs ##

M2_pred <- summary(M2)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M2_pred_ci09 <- summary(M2_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_2 <- as.data.frame(M2_pred)
data_figure_2_09 <- as.data.frame(M2_pred_ci09)

data_figure_2 <- cbind(data_figure_2, data_figure_2_09[,3:4])

data_figure_2 <- data_figure_2  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_2 <- rbind(data_figure_2, M2_placebo_data)

data_figure_2[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1b <- ggplot(data = data_figure_2, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Autocratic Consolidation", 
       y = "Estimated effect of Autocratic Consolidation", 
       title = "Mahalanobis: Terrorist Attacks")


#### PS Match #### 

PM_M3 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "regressed_autocracy", refinement.method = "ps.match", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M3 <- PanelEstimate(sets = PM_M3, data = vdem)

M3_ci09 <- PanelEstimate(sets = PM_M3, data = vdem, confidence.level = .9)

M3_placebo <- placebo_test(PM_M3,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M3_placebo_t_2_est <- M3_placebo$estimates[1]
M3_placebo_t_2_high <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.975)
M3_placebo_t_2_low <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.025)

M3_placebo_t_2_high_90 <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.95)
M3_placebo_t_2_low_90 <- quantile(M3_placebo$bootstrapped.estimates[ ,1],0.05)


M3_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M3_placebo_t_2_est,
                              conf_low_95 = M3_placebo_t_2_low,
                              conf_low_90 = M3_placebo_t_2_low_90,
                              conf_high_95 = M3_placebo_t_2_high, 
                              conf_high_90 = M3_placebo_t_2_high_90)

#95% CIs ##

M3_pred <- summary(M3)

M3_pred <- M3_pred$summary

M3_pred <- cbind(M3_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M3_pred_ci09 <- summary(M3_ci09)

M3_pred_ci09 <- M3_pred_ci09$summary

M3_pred_ci09 <- cbind(M3_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_3 <- as.data.frame(M3_pred)
data_figure_3_09 <- as.data.frame(M3_pred_ci09)

data_figure_3 <- cbind(data_figure_3, data_figure_3_09[,3:4])

data_figure_3 <- data_figure_3  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_3 <- rbind(data_figure_3, M3_placebo_data)

data_figure_3[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1c <- ggplot(data = data_figure_3, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Autocratic Consolidation", 
       y = "Estimated effect of Autocratic Consolidation", 
       title = "PS Match: Terrorist Attacks")




#### CBPS Match #### 

PM_M4 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "regressed_autocracy", refinement.method = "CBPS.match", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M4 <- PanelEstimate(sets = PM_M4, data = vdem)

M4_ci09 <- PanelEstimate(sets = PM_M4, data = vdem, confidence.level = .9)

M4_placebo <- placebo_test(PM_M4,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M4_placebo_t_2_est <- M4_placebo$estimates[1]
M4_placebo_t_2_high <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.975)
M4_placebo_t_2_low <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.025)

M4_placebo_t_2_high_90 <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.95)
M4_placebo_t_2_low_90 <- quantile(M4_placebo$bootstrapped.estimates[ ,1],0.05)


M4_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M4_placebo_t_2_est,
                              conf_low_95 = M4_placebo_t_2_low,
                              conf_low_90 = M4_placebo_t_2_low_90,
                              conf_high_95 = M4_placebo_t_2_high, 
                              conf_high_90 = M4_placebo_t_2_high_90)

#95% CIs ##

M4_pred <- summary(M4)

M4_pred <- M4_pred$summary

M4_pred <- cbind(M4_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M4_pred_ci09 <- summary(M4_ci09)

M4_pred_ci09 <- M4_pred_ci09$summary

M4_pred_ci09 <- cbind(M4_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_4 <- as.data.frame(M4_pred)
data_figure_4_09 <- as.data.frame(M4_pred_ci09)

data_figure_4 <- cbind(data_figure_4, data_figure_4_09[,3:4])

data_figure_4 <- data_figure_4  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_4 <- rbind(data_figure_4, M4_placebo_data)

data_figure_4[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1d <- ggplot(data = data_figure_4, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Autocratic Consolidation", 
       y = "Estimated effect of Autocratic Consolidation", 
       title = " CBPS Match: Terrorist Attacks")



#### CBPS Weight #### 

PM_M5 <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                    treatment = "regressed_autocracy", refinement.method = "CBPS.weight", 
                    data = vdem, match.missing = FALSE, 
                    covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                      I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                      I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                      I(lag(e_area_log, lag_2)) + 
                      I(lag(ma_count_attacks_dom, lag_2)), 
                    size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                    lead = lead_5, forbid.treatment.reversal = FALSE,
                    use.diagonal.variance.matrix = TRUE, placebo.test = T)

M5<- PanelEstimate(sets = PM_M5, data = vdem)

M5_ci09 <- PanelEstimate(sets = PM_M5, data = vdem, confidence.level = .9)

M5_placebo <- placebo_test(PM_M5,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M5_placebo_t_2_est <- M5_placebo$estimates[1]
M5_placebo_t_2_high <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.975)
M5_placebo_t_2_low <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.025)

M5_placebo_t_2_high_90 <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.95)
M5_placebo_t_2_low_90 <- quantile(M5_placebo$bootstrapped.estimates[ ,1],0.05)


M5_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M5_placebo_t_2_est,
                              conf_low_95 = M5_placebo_t_2_low,
                              conf_low_90 = M5_placebo_t_2_low_90,
                              conf_high_95 = M5_placebo_t_2_high, 
                              conf_high_90 = M5_placebo_t_2_high_90)

#95% CIs ##

M5_pred <- summary(M5)

M5_pred <- M5_pred$summary

M5_pred <- cbind(M5_pred , time = c(0,1,2,3,4,5))

# 90% CIs ##

M5_pred_ci09 <- summary(M5_ci09)

M5_pred_ci09 <- M5_pred_ci09$summary

M5_pred_ci09 <- cbind(M5_pred_ci09 , time = c(0,1,2,3,4,5))

data_figure_5 <- as.data.frame(M5_pred)
data_figure_5_09 <- as.data.frame(M5_pred_ci09)

data_figure_5 <- cbind(data_figure_5, data_figure_5_09[,3:4])

data_figure_5 <- data_figure_5  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_5 <- rbind(data_figure_5, M5_placebo_data)

data_figure_5[nrow(data_figure_1) + 1,] <- c(-1, 0, NA,  NA, NA, NA)

#### Building Subplots for each Estimand ####

f1e <- ggplot(data = data_figure_5, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Autocratic Consolidation", 
       y = "Estimated effect of Autocratic Consolidation", 
       title = "CBPS Weight: Terrorist Attacks")

##### GG Arrange Plots #####

ggarrange(f1a, f1b, f1c, f1d, f1e, nrow = 5)

ggsave("outputs/Figure_G5.png", units= c("cm"), width = 20, height = 35, dpi = 900)
ggsave("outputs/Figure_G6.pdf", units= c("cm"), width = 20, height = 35, dpi = 900)

